# Paper: Is it still the economy? Economic voting in polarized politics
# Individual-level simulation
# Replication: Figure S8 (Supplementary Material)
remove(list = ls())

library(conflicted)
library(tidyverse)
library(dplyr)
conflict_prefer("filter","dplyr")
library(ggplot2)
library(here) # When we have an R-project, it reads files without setting an directory
conflict_prefer("here", "here")
library(nnet)
library(ggeffects)
library(MASS)
conflict_prefer("select","dplyr")

load(here("Final_Paper", "R&R", "IndLevel.RData"))


rep_datar <- model_data %>% filter(rep_inc == 1) %>% 
  filter(ptd_id3 == "Rep")
table(rep_datar$state, useNA = "always")

summary(ml_repr <- multinom(vote_inc_choice ~ econ_eval*percep_polar + 
                              female + black + hispanic + asian + other_race + 
                              log_age + college + income + state_fips + 
                              election, 
                            data = rep_datar))
table(model_data$percep_polar)
ml_repr2 <- ml_repr %>% 
  ggeffect(terms = c("econ_eval", "percep_polar[0.1, 0.5, 1]"), 
           ci.lvl = 0.95)
glimpse(ml_repr2)
ml_repr2 <- data.frame(ml_repr2)

ml_repr2 <- ml_repr2 %>% 
  mutate(ptd_id = "Rep") %>% 
  mutate(incumbent = "Rep")

rep_datad <- model_data %>% filter(rep_inc == 0) %>% 
  filter(ptd_id3 == "Rep")


summary(ml_repd <- multinom(vote_inc_choice ~ econ_eval*percep_polar + 
                              female + black + hispanic + asian + other_race + 
                              log_age + college + income + state_fips + 
                              election, 
                            data = rep_datad))
table(model_data$percep_polar)
ml_repd2 <- ml_repd %>% 
  ggeffect(terms = c("econ_eval", "percep_polar[0.1, 0.5, 1]"), 
           ci.lvl = 0.95)
glimpse(ml_repd2)
ml_repd2 <- data.frame(ml_repd2)

ml_repd2 <- ml_repd2 %>% 
  mutate(ptd_id = "Rep") %>% 
  mutate(incumbent = "Dem")

ml_rep <- rbind(ml_repr2, ml_repd2)
glimpse(ml_rep)

g_rep_bad <- ml_rep %>% 
  mutate(response.level = factor(response.level,
                                 levels = c("Inc", 
                                            "Ops",
                                            "Abs",
                                            "other"),
                                 labels = c("Incumbent", 
                                            "Oppostion",
                                            "Abstention",
                                            "Other"))) %>% 
  mutate(incumbent = factor(incumbent,
                            levels = c("Rep", "Dem"))) %>% 
  mutate(group = as.factor(group)) %>% 
  #mutate(response.level = as.factor(response.level)) %>%
  filter(x == 1) %>% 
  mutate(x = as.factor(x)) %>%
  #filter(response.level == "Inc") %>% 
  #filter(group == 0.33 | group == 1) %>% 
  ggplot(aes(colour = incumbent, shape = group, fill = incumbent)) + 
  #geom_hline(yintercept = 0, colour = "black", lty = 2, size = .75) +
  geom_pointrange(aes(x = response.level, y = predicted, 
                      ymin = conf.low, ymax = conf.high),
                  lwd = .65, position = position_dodge(width = .5)) + 
  theme_light() + 
  labs(x = NULL, y = "Predicted Probability", 
       title = "(a) Republicans") +
  scale_color_manual(name = "Presidency:", 
                     values = c("Red",  "Blue"),
                     labels = c("Rep", "Dem")) +
  scale_shape_manual(name = "Ideological gap:",
                     labels = c("0.1", ".5", "1"),
                     values = c(22, 21, 24)) +
  scale_fill_manual(name = "Presidency:",
                    values = c("Red",  "Blue"), 
                    labels = c("Rep", "Dem")) +
  scale_y_continuous(breaks = c(1, 0.8, 0.6, 0.4, 0.2, 0),
                     limits = c(0,1)) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 12),
        axis.text.y = element_text(size = 10),
        axis.title.y = element_text(size = 14),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 11),
        legend.position = "none",
        legend.direction = "horizontal",
        plot.title = element_text(size = 15))
g_rep_bad


####################################################
####################################################
####################################################
######### Democratic Voters
####################################################
####################################################
####################################################

dem_datar <- model_data %>% filter(rep_inc == 1) %>% 
  filter(ptd_id3 == "Dem")

summary(ml_demr <- multinom(vote_inc_choice ~ econ_eval*percep_polar + 
                              female + black + hispanic + asian + other_race + 
                              log_age + college + income + state_fips +
                              election, 
                            data = dem_datar))
table(model_data$percep_polar)
ml_demr2 <- ml_demr %>% 
  ggeffect(terms = c("econ_eval", "percep_polar[0.1, 0.5, 1]"), 
           ci.lvl = 0.95)
glimpse(ml_demr2)
ml_demr2 <- data.frame(ml_demr2)

ml_demr2 <- ml_demr2 %>% 
  mutate(ptd_id = "Dem") %>% 
  mutate(incumbent = "Rep")

dem_datad <- model_data %>% filter(rep_inc == 0) %>% 
  filter(ptd_id3 == "Dem")

table(dem_datad$econ_eval)
summary(ml_demd <- multinom(vote_inc_choice ~ econ_eval*percep_polar + 
                              female + black + hispanic + asian + other_race + 
                              log_age + college + income + state_fips + 
                              election, 
                            data = dem_datad))
table(model_data$percep_polar)
ml_demd2 <- ml_demd %>% 
  ggeffect(terms = c("econ_eval", "percep_polar[0.1, 0.5, 1]"), 
           ci.lvl = 0.95)
glimpse(ml_demd2)
ml_demd2 <- data.frame(ml_demd2)

ml_demd2 <- ml_demd2 %>% 
  mutate(ptd_id = "Dem") %>% 
  mutate(incumbent = "Dem")

ml_dem <- rbind(ml_demr2, ml_demd2)
glimpse(ml_dem)

g_dem_bad <- ml_dem %>% 
  mutate(response.level = factor(response.level,
                                 levels = c("Inc", 
                                            "Ops",
                                            "Abs",
                                            "other"),
                                 labels = c("Incumbent", 
                                            "Oppostion",
                                            "Abstention",
                                            "Other"))) %>% 
  mutate(incumbent = factor(incumbent,
                            levels = c("Rep", "Dem"))) %>% 
  mutate(group = as.factor(group)) %>% 
  #mutate(response.level = as.factor(response.level)) %>%
  filter(x == 1) %>% 
  mutate(x = as.factor(x)) %>%
  #filter(response.level == "Inc") %>% 
  #filter(group == 0.33 | group == 1) %>% 
  ggplot(aes(colour = incumbent, shape = group, fill = incumbent)) + 
  #geom_hline(yintercept = 0, colour = "black", lty = 2, size = .75) +
  geom_pointrange(aes(x = response.level, y = predicted, 
                      ymin = conf.low, ymax = conf.high),
                  lwd = .65, position = position_dodge(width = .5)) + 
  theme_light() + 
  labs(x = NULL, y = "Predicted Probability", 
       title = "(b) Democrats") +
  scale_color_manual(name = "Presidency:", 
                     values = c("Red",  "Blue"),
                     labels = c("Rep", "Dem")) +
  scale_shape_manual(name = "Ideological gap:",
                     labels = c("0.1", ".5", "1"),
                     values = c(22, 21, 24)) +
  scale_fill_manual(name = "Presidency:",
                    values = c("Red",  "Blue"), 
                    labels = c("Rep", "Dem")) +
  scale_y_continuous(breaks = c(1, 0.8, 0.6, 0.4, 0.2, 0),
                     limits = c(0,1)) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 12),
        axis.text.y = element_text(size = 10),
        axis.title.y = element_text(size = 14),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 11),
        legend.position = "none",
        legend.direction = "horizontal",
        plot.title = element_text(size = 15))
g_dem_bad


####################################################
####################################################
####################################################
######### Independents
####################################################
####################################################
####################################################

ind_datar <- model_data %>% filter(rep_inc == 1) %>% 
  filter(ptd_id3 == "Ind")

summary(ml_indr <- multinom(vote_inc_choice ~ econ_eval*percep_polar + 
                              female + black + hispanic + asian + other_race + 
                              log_age + college + income + state_fips + 
                              election, 
                            data = ind_datar))
table(model_data$percep_polar)
ml_indr2 <- ml_indr %>% 
  ggeffect(terms = c("econ_eval", "percep_polar[0.1, 0.5, 1]"), 
           ci.lvl = 0.95)
glimpse(ml_indr2)
ml_indr2 <- data.frame(ml_indr2)

ml_indr2 <- ml_indr2 %>% 
  mutate(ptd_id = "Ind") %>% 
  mutate(incumbent = "Rep")

ind_datad <- model_data %>% filter(rep_inc == 0) %>% 
  filter(ptd_id3 == "Ind")


summary(ml_indd <- multinom(vote_inc_choice ~ econ_eval*percep_polar + 
                              female + black + hispanic + asian + other_race + 
                              log_age + college + income + state_fips + 
                              election, 
                            data = ind_datad))
table(model_data$percep_polar)
ml_indd2 <- ml_indd %>% 
  ggeffect(terms = c("econ_eval", "percep_polar[0.1, 0.5, 1]"), 
           ci.lvl = 0.95)
glimpse(ml_indd2)
ml_indd2 <- data.frame(ml_indd2)

ml_indd2 <- ml_indd2 %>% 
  mutate(ptd_id = "Ind") %>% 
  mutate(incumbent = "Dem")

ml_ind <- rbind(ml_indr2, ml_indd2)
glimpse(ml_ind)

g_ind_bad <- ml_ind %>% 
  mutate(response.level = factor(response.level,
                                 levels = c("Inc", 
                                            "Ops",
                                            "Abs",
                                            "other"),
                                 labels = c("Incumbent", 
                                            "Oppostion",
                                            "Abstention",
                                            "Other"))) %>% 
  mutate(incumbent = factor(incumbent,
                            levels = c("Rep", "Dem"))) %>% 
  mutate(group = as.factor(group)) %>% 
  #mutate(response.level = as.factor(response.level)) %>%
  filter(x == 1) %>% 
  mutate(x = as.factor(x)) %>%
  #filter(response.level == "Inc") %>% 
  #filter(group == 0.33 | group == 1) %>% 
  ggplot(aes(colour = incumbent, shape = group, fill = incumbent)) + 
  #geom_hline(yintercept = 0, colour = "black", lty = 2, size = .75) +
  geom_pointrange(aes(x = response.level, y = predicted, 
                      ymin = conf.low, ymax = conf.high),
                  lwd = .65, position = position_dodge(width = .5)) + 
  theme_light() + 
  labs(x = NULL, y = "Predicted Probability", 
       title = "(c) Independents") +
  scale_color_manual(name = "Presidency:", 
                     values = c("Red",  "Blue"),
                     labels = c("Rep.", "Dem.")) +
  scale_shape_manual(name = "Ideological gap:",
                     labels = c("0.1", "0.5", "1"),
                     values = c(22, 21, 24)) +
  scale_fill_manual(name = "Presidency:",
                    values = c("Red",  "Blue"), 
                    labels = c("Rep.", "Dem.")) +
  scale_y_continuous(breaks = c(1, 0.8, 0.6, 0.4, 0.2, 0),
                     limits = c(0,1)) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 12),
        axis.text.y = element_text(size = 10),
        axis.title.y = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        legend.position = "bottom",
        legend.direction = "horizontal",
        plot.title = element_text(size = 15))
g_ind_bad


multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  if (numPlots==1) {
    print(plots[[1]])
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}
multiplot(g_rep_bad, g_ind_bad, g_dem_bad, cols = 2)



